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ABSTRACT 



A computer simulation model is developed treating wave 
propagation in a random, inhomogeneous ocean as transmission 
through a linear time-invariant, space-variant random communi- 
cation channel. The ocean volume is modelled by an index of 
refraction which is decomposed into a depth-dependent deter- 
ministic part and a depth-independent Gaussian zero-mean 
random part. Computer simulated output electrical signals 
were generated that depend on the complex frequency spectrum 
of the transmitted electrical signal, the far-field beam 
pattern of the transmit array and the random transfer function 
of the ocean medium. Output was generated for different test 
cases. In all cases the transmit electrical signal was repre- 
sented by a finite Fourier series and random cases were 
modelled by a random number generator. The computer simulated 
output electrical signals were then processed by a 3-D DFT 
beamformer and the results for the deterministic inhomogeneous 
and random inhomogeneous cases were compared to the homogeneous 
non-random case in order to study the effects of the medium 
on signal distortion and source localization. 
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INTRODUCTION 



I . 

Based on the linearized acoustic wave equation for small 
amplitude acoustic signals, sound propagation through the ocean 
medium can be modelled as a signal passing through a linear 
time-variant , space-variant random filter. The term "space- 
variant" implies that the sound speed profile is a function 
of position. The space variant property results in scatter 
or angular spread due to refraction. If the filter is 
"space-variant," then an isospeed medium is implied. As a 
result, there will be no refraction, and hence, no scatter or 
angular spread since the sound rays will be travelling in 
straight lines. The term "time-variant" implies motion of, 
for example, discrete point scatterers , the ocean surface, 
and the transmit and receive arrays (apertures) . The time- 
variant property results in both Doppler spread and spread 
in time delay. 

By using a linear system's theory approach, different 
propagation paths, such as direct paths, bottom reflected 
paths and surface reflected paths, can be modelled as the 
outputs from linear filters. Furthermore, different trans- 
mit signals and transmit and receive far-field directivity 
functions can easily be coupled to various transfer functions 
of the random, inhomogeneous ocean medium in order to study 
their effects on signal distortion and source localization 
using various space- time signal processing algorithms. 
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Work based on this approach has been done by Laval [Refs. 
5,6], Laval and Labasque [Ref. 7], Middleton [Refs. 8,9] 
and Ziomek [Ref. 1]. Middleton [Refs. 8,9] described propaga- 
tion phenomena in a random inhomogeneous ocean with the use 
of space-time operators but did not concern himself directly 
with the derivation of random, time-variant, space-variant 
transfer functions. 

Ziomek [Ref. 2] derived a time-invariant space-variant 
random transfer function of the ocean volume based on the WKB 
approximation. Ziomek [Ref. 3] also derived an equation for 
the random, output electrical signal at each element in a 
receive planar array of complex weighted point sources in 
terms of the frequency spectrum of the transmitted electrical 
signal, the transmit and receive arrays, and the random ocean 
medium transfer function. 

The purpose of this thesis is to implement on a computer 
the equation for the output electrical signal as derived by 
Ziomek [Ref. 3]. This mathematical computer model simulates 
the effects of wave propagation through a random inhomogeneous 
ocean. The ocean volume is modelled by the transfer function 
derived in Ziomek [Ref. 2] and incorporates an index of 
refraction which is a function of depth. The index of refrac- 
tion is decomposed into a depth-dependent deterministic com- 
ponent and a depth-independent zero mean Gaussian random 
component . 

Section II is devoted to a theoretical description of the 
computer model. The notation and geometry of the computer 
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model are explained and the fundamental equations on which 
the computer model is based, are stated. Also, some additional 
constraints are made to allow for simpler implementation. 

Section III describes the numerical techniques used and 
the resulting mathematical procedures which are implemented 
in the program. Furthermore, a brief outline of the computer 
program is given. 

In Section IV a number of test cases are defined which 
were used to test the computer model. Computer simulated 
output electrical signals were generated for test cases 
incorporating deterministic homogeneous, deterministic 
inhomogeneous and random inhomogeneous medium models. The 
output signals were processed by a 3-D DFT beamformer and 
the results for the deterministic homogeneous case were used 
as the baseline case. Both deterministic and random 
inhomogeneous cases were compared to this baseline case in 
order to study the effects of the medium on signal distortion 
and source localization. 
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II. THEORY ON THE COMPUTER SIMULATION MODEL 



A. MODEL DESCRIPTION 

The mathematical computer model calculates the output 
electrical signal at each element of a receive planar array 
of point sources in terms of the complex frequency spectrum 
of the transmitted electrical signal, the transmit and receive 
arrays, and the random ocean medium transfer function. The 
theory for the model is based on time-variant, space-variant 
random filters as discussed in Ziomek [Ref. 1]. Ziomek [Ref. 

2] derived a time-invariant space-variant transfer function 
for a random ocean volume where the index of refraction, which 
is composed of a deterministic and a random component, is a 
function of depth. This transfer function is based on the 
WKB approximation, which is a solution to the wave equation 
when the index of refraction is a function of depth. Further- 
more, Ziomek [Ref. 3] derived an equation for the random 
output electrical signals appearing at each element in a 
receive planar array of complex weighted point sources in 
terms of the frequency spectrum of the transmitted electrical 
signal, the transmit and receive arrays, and the random ocean 
transfer function. The computer model combines and implements 
the equations derived by Ziomek [Refs. 2,3] for the transfer 
function and the output electrical signals. 

1 . Geometry of the Model 

Figure 1 shows the geometry of the communication 
channel as simulated by the computer model. The model 
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Fig. 1 Model Geometry 



describes the wave propagation between a transmit planar array 
and a receive planar array whose centers are located at 



(x o ,y o ,z Q ) and ( x r 'y r < z r ) ' respectively. The x, y and z axes 
describe cross range, depth and range, respectively. The 
transmit array, which is a parallel to the XY plane, consists 
of M xN complex weighted point sources, where M and N are odd. 
The sources are equally spaced in the X- and Y-directions 
with spacings of d^ and d^ meters, respectively. Note that 
in general d^ ^ d^ . The complex weights are assumed to be 
separable and given by 



c = a exp { j 0 } 
m m c J m 



( 2 . 1 ) 



for the weighting in the X-direction and by 



d n = b n exp (2 • 2) 

for the weighting in the Y-direction. The combined weight 

factor for each transmit element is c d , where m and n are 

m n 

the indices describing the transmit element at position 

(x +md ,y +nd ,z ). The complex weights are used for ampli- 
o x o y o 

tude shading and beamsteering of the transmit beam pattern. 

The receive planar array, which is also parallel to the XY 

plane, consists of M’ *N' elements, where M' and N' are odd. 

The sources are equally spaced in the X- and Y-directions with 

spacings of d' and d' meters, respectively, where in general 
x y 

d' ^ d'. The position of a receiver element with indices 
x y 

(i,q) is given by (x +id ' ,y +qd ’ , z ) . 

J. 2 S. J_ y L 
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The speed of sound in the medium is assumed to be 
only a function of depth and is given by c(y) . The reference 
speed of sound c q is the speed of sound at the center of the 
transmit array, i.e., c q = c(y Q ). The index of refraction 
n(y) is a function of depth and is composed of a determinis- 
tic component and a random component: 



n(y) = C Q /c(y) = n D (y) + n R (y) ( 2 . 3 ) 

where n is the deterministic component and n_ is the random 
component of the index of refraction. The initial propagation 
vector will be denoted by k n , k Q = 2Trf/c o is the magnitude 
of the initial propagation vector at depth y = y Q (i.e., at 
the center of the transmit array) . The components of the 
initial propagation vector along the x, y and z axes are 
described in terms of the direction cosines u , v and w : 



k = k u 
x o o 



k = k v 
y o o 



k = k w 
z o o 



and 



2 .22 

w = 1 - u - V 

o o o 



( 2 . 4 ) 

( 2 . 5 ) 

( 2 . 6 ) 

( 2 . 7 ) 



Note that the magnitude of the initial propagation vector at 
the transmit elements which are located at a depth nd meters 

y 
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from the center of the transmit array is, in general, not 
equal to k Q and is given by : 

k n = 27rf/c(y o +nd y ) (2.8) 

One can also use spatial frequencies (in cycles/meter) to 
denote directions of propagation: 

f = u /X; f = v /X; f = w /X (2.9) 

x o y o' z o' 



2 . Transfer Function 

The medium transfer function is given by a function of 
frequency f, spatial frequency f and receiver depth y for 
a given source depth y Q [Ref. 2] : 



H M (f ' f y ;Y) = A M (f,f y ;Y)eXp{j [0 MD (f ' f y ;y) 



+ 9 MR (f 'V y)]} 



( 2 . 10 ) 



where is the amplitude term given by 



A M - <2*V' 1/2 = <k o V o rV2 (2 - 11: 



and 0 are the deterministic and random phase terms, 

respectively : 



MD 



(-k /4irf ) 

o y 



[n D (?) - l]d? 



( 2 . 12 ) 



' 16 



(2.13) 



MR 



2 , Y 

(-k 0 /2uf ) / n D U)n R U)d? 



This transfer function is valid in the absence of turning 
points and when the following inequality is obeyed: 



n 2 (y) - 1 | / v 2 < 1 



(2.141 



In the case of a speed of sound profile with a constant 
gradient g, the deterministic component of the index of 
refraction becomes: 



n n(y) = c n/ + g(y-yj ) 



(2.151 



Substitution into equation (2.12) results in 



MD 



k C 

2^T { V tn D (Y) _1] + y " y o } 



(2. 16] 



As far as the random phase term is concerned we will 

proceed by assuming that the random component of the index of 

refraction is not a function of depth and can be represented 

2 

by a Gaussian, zero mean random variable with variance a : 



n R ( C ) = n R = an NR 



(2.171 



where n is N(0,a ) and n which is the normalized random 
K NK 

component of the index of refraction, is N(0,1) . When n R 
is not a function of depth, equation (2.13) becomes: 
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(2.18) 



9 MR 



k 

o 



V 

o 



n 



R 



/ 



y 



y 



o 



n D (C)de 



Substituting the equation of n D for a constant gradient into 
equation (2.18) leads to: 



0 MR 



V g n R ln ^ n D * 
o ^ 



(2.19) 



For a homogeneous medium, equations (2.16) and (2.19) become: 



0 MD 



0 



( 2 . 20 ) 



9 MR - -^ n R (y - y o } (2 ' 21) 

o 

Equations (2.11) through (2.13) for the amplitude and 
phase terms of the transfer function show that the medium has 
the effect of angle modulating the transmitted sound field, 
also referred to as scattering. 

3 . Output Electrical Signal 

The output electrical signal at a receive location 
(x r +id^,y r +qd^,z r ) is given by [Ref. 3] : 



oo 11 

y(t,x +id^,y +qd',z ) = / X(f) / / f 2 • 

* -co -la 

q 



(N-l) /2 

l 

n=- (N-l) /2 
(M-l)/2 

l 

m=-(M _ l)/2 



( V c X (f 'V Y) ^ { 'Vo% }e ^ { ' jk n (1 ‘V v o )1/2AZl * 

c exp{-jk u AX. }dv du exp{ j2TTft}df ; u 2 +v 2 < 1 (2.22) 

m norm o o ^ J oo — 
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where 



a q = Un2(y r+ qd^) -1|] 1/2 



AY = y - y + qd ' - nd 
qn J r J o ^ y y 



AX. = x - x + id' - nd 
1m r o x x 



AZ = z - z 
r o 



X(f) : complex frequency spectrum of the transmitted 

electrical signal. 



c n = c (y Q +n dy ) ' speed of sound at transmit 
element with index n. 



k = 27 rf/c(y + nd v ) ; magnitude of initial propagation 
vector at transmit element with index n. 



m 



n 



complex weights of the transmit array. 



H (f,f ;y) : Random time-invariant space-variant 

^transfer function of the ocean medium. 



Note that H„. is a function of both indices 
q and n, i.e., it is a function of the y 
coordinate of the source location y Q + nd 
and receiver location y + qd ' . This can^ 
also be expressed by using the notation 



H (k n' v 0 ' n,q) 



wave number k 
indices n 



expressing H as a function of 



and n q . 



M 

direction cosine v and 



The above equations show how the medium transfer 
function is used in describing the acoustic field at a receiver 
location as a function of the acoustic field at a transmit 
location. The far-field directivity functions of transmit and 



19 



receive arrays are used to couple an electrical signal to 
the medium and vice-versa. 

The region of integration over direction cosine v q 
is limited by a^ . This expresses the fact that the transfer 
function is not valid for direction cosine values v q approach- 
ing 0 (close to a turning point) . In most applications the 
limits of integration for both integrals can be further limited 
as shown in Section III. 

B. CONSTRAINTS ON THE MODEL 

In this section the inherent constraints on the model will 
be investigated. Also, some additional constraints will be made 
1 . Geometrical Constraints 

These constraints concern the transfer function which 
is only valid under the two following restrictions : 

a. The assumption made in the derivation of the transfer 
function, in order to validate a binomial expansion was [Ref. 2] 

| [n 2 (y) -l]/v 2 ] | < 1 (2.23) 

Thus before using the computer model one should evalu- 
ate equation (2.23) for the chosen sound speed profile and 
array locations to make sure the inequality is obeyed. 

b. Absence of turning points. The occurrence of turning 
points in a particular transmitter/receiver geometry depends 
on the initial direction of propagation, described by the 
angle 6 q between the initial propagation vector and the 
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Y-axis. Using Snell's law we find that at the turning point 



depth y T : 



c(y T } = C o /sin(e o ) 



(2.24) 



The relation between the angle and v q is given by: 



v 



o 



cos (6 ) 
o 



(2.25) 



Combining equations (2.24) and (2.25) into an inequality 
which has to be obeyed to avoid a turning point results in: 



a sound speed profile which is assumed to have a constant 
gradient g and a random component of the index of refraction 
n^ which is not a function of depth. Note that the equations 
derived in Ziomek [Refs. 2,3] are more general and that the 
computer program described in this thesis can be easily modi- 
fied to handle a more complex sound speed profile. For purpose 
of computer simulation, the random component of the index of 

refraction was assumed to be a zero mean, Gaussian random 

2 

number with variance a and was modelled by a random number 




(2.26) 



and we conclude that obeying the inequality (2.23) also 
guarantees the absence of turning points. 



2 . Sound Speed Profile 

The results presented in this thesis are based upon 
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generator. These restrictions make it possible to evaluate 
the transfer function with the closed form expressions given 
in equations (2.11), (2.16) and (2.19) for amplitude, deter- 

ministic phase component and random phase component of the 
transfer function, respectively. In the calculation of the 
random phase terms of the transfer function, a different 
random number was drawn for each combination of transmit 
element depth y Q + nd^ and receive element depth y^_ + qd^ , 
thereby assuming that the random phase terms for the different 
transmission paths are uncorrelated. 

3 . Frequency Spectrum of the Transmitted Electrical 
Signal 

The transmitted electrical signal was represented by 

a finite Fourier series with fundamental period T q and maximum 

frequency f Hz: 

^ max 

f = KMAX • (2.27) 

max T 

o 

where KMAX denotes the total number of harmonics in the signal. 
This representation does not impose a severe restriction on 
the input signal, since every finite energy signal can be 
represented in the sense of a least mean-square error by a 
finite Fourier series. The expression for the frequency 
spectrum becomes : 



KMAX 

X(f) = l [c, 6 ( f - kf ) + c * 6 ( f + k f )] (2.28) 

-i _ K O K O 
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where f is the fundamental frequency 1 /T q and the DC term 
is assumed to be zero. The complex Fourier series coefficients 
c^ determine the amplitude and phase relations among the 
different frequency components. The integral in equation 
(2.22) now reduces to a summation and the expression for the 
real output electrical signal at a receiver element with 
indices i and q becomes : 



y (t,i,q) = 2 Re { l c,(kf ) 2 / / • 

k=l -1 a 

q 



KMAX _ 1 1 



KMAX 




(N-l)/2 



(M-l)/2 



m=-(M-l)/2 



c exp [-jk u AX. ]dv du exp[j2TTkf t]}; 
m^ J normoo^ J o 



(2.29) 
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Ill . 



IMPLEMENTATION OF THE MODEL 



The main task of the computer simulation model is to 
compute and file the electrical output signal at all elements 
of the receive planar array according to equation (2.29) . The 
model was implemented on an IBM system/370 mainframe computer 
using FORTRAN. The program is computation intensive because 
the evaluation of equation (2.29) involves a double integral 
with nested summations. Therefore, much attention has been 
given to program efficiency. This section discusses the 
equations which were actually implemented and the overall logic 
of the program. 

A. NUMERICAL METHODS 

The model computes the real output electrical signal at 
all receiver elements as given by equation (2.29). Note that 
the double integral with respect to (w.r.t.) direction 
cosines u q and v q behaves like a frequency response or trans- 
fer function H(f,i,q) , which depends on the position of the 
receiver elements as indicated by the indices (i.q) . We define 



A f 



.2 



1 1 



(N-l)/2 



H(f ,i,q) 




/ / 




a 



q 



exp[-jk v AY ]exp[-jk (1-u^-v^) 
* J n o qn ^ J n o o 



(M-l)/2 



no qn 



'AZ] l 

m=-(M-l)/2 



c exp[-jk AX. u ]dv du 
m^ J nimo oo 



(3.1) 
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2 

Note that the factor l/ c n i n equation (2.29) was moved in 

2 

front and replaced by the factor 1/c q . This is possible 
because the differences in sound speed at the different trans- 
mitter elements are very small for all practical cases and 
the effect on the amplitude of H(f,i,q) can be neglected. 
However, the dependence of the initial speed of sound on the 
index n of the transmit array will be taken into account in 
the evaluation of the phase terms as expressed by the different 
values of k n in equation (3.1) . 

Using the definition of H(f,i,q), the expression for the 
real output electrical signal at a receiver element with 
indices (i,q) becomes 

KMAX 

y(t,i,q) = 2 Re { l c.H (kf , i,q) exp ( j 2irkf t) } (3.2) 

k=l K ° ° 

The computer program model first calculates H(f,i,q) 
as given in equation (3.1). Then, with the use of H(f,i,q), 
the output electrical signal is readily computed according 
to equation (3.2). 

The ocean medium transfer function is possibly random 
and therefore difficult to integrate numerically so a change 
in the order of integration in equation (3.1) is made to 
achieve a more suitable form for implementation: 
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H(f ,i,q) 



2 1 (N-l)/2 

= — / 7 d H..(k ,v ,n,q)exp[-jk v AY ] 

2 1 /0 n M n o M ^ J n o qn 

c a n=-(N-l)/2 ^ 

o q 



1 (M-l) /2 _ 9 i /p 

/ £ c exp[-jk u AX. ] exp [— jlc (l-u^-vj 7 AZ]du dv (3.3) 

J m n o un n oo r\ 



-1 m=-(M-l )/2 



o o 



In order to improve efficiency, we concentrate first on 
the inner integral w.r.t. u q , which we rewrite as 



-1 



/ S v ( f , u ) exp{ -k [u AX. + (l-u 2 -v 2 ) 1 // 2 AZ] }du (3.4) 

'Xo noi oo o 



where 



AX . 

l 



- x 



+ 



id* 

x 



Note that 



S X (f ' u o } = 



(M-D/2 

l 

m =- ( m - 1)/2 



c exp[jk u md ] 
m ^ J n o x 



(3.5) 



is the transmit far-field directivity function of a line array 
consisting of equally spaced point-sources as a function of 
direction cosine u q and frequency f. Now equation (3.3) reduces 
to : 

2 1 (N-l)/2 

H(f ,i,q) - -j / I W k n' v o' n ' q)ex P 1 ' jk n v o A V ' 

C Q a n=-(N-l)/2 ^ 

/ S (f,u )exp{-jk [u AX. + (l-u 2 -v 2 ) ^ 2 AZ] }du dv (3.6) 

, J x' ' o ^ J n o l oo oo 
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In the case of simple amplitude weighting in the x-direction 
such as rectangular, triangular or Hamming weighting, we can 
obtain closed form expressions for S (f,u ) . This improves 
program efficiency. For example, in the case of rectangular 
amplitude weighting 



sm [k (u - u, )Md /2] 

S (f u ) = - ° b X 

^X ir ' o' sin [k (u -u,)d /2] 

n o b x' 



(3.7) 



where u, is the direction cosine value to which the transmit 
b 

far-field beam pattern is steered. 

Before presenting more details on how H(f,i,q) is calcu- 
lated, a brief overview on numerical and analytical integration 
techniques which were used is given. 

1 . Integration Techniques 

a. Method of Stationary Phase 

This method approximates an integral of the form 

b 

I / S(z) exp[jg(z)]dz (3.8) 

a 



where g(z) is a rapidly varying phase term compared to the 
slowly varying amplitude term S(z) . Officer [Ref. 10] shows 
that this integral can be approximated by using the fact that 
the main contribution exists at the stationary point z g ^ which 
is defined by 



ag<V 

— ap- = 0 



(3.9) 
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The limits of the integral may be extended to infinity 
because only the region around the stationary point contributes 
to the integral. Assume that around the stationary point, z 
is given by 



z = z + £ (3.10) 

sp 

The approximation of the integral is then given by 

1 = S(z sp ) ^ 27T / I g" (z sp ) I exp{j [g (z gp ) ±J] } (3.11) 

under the condition: 



kg'" (z ) | 



<< l 



(3.12) 



where the sign in the phase term of equation (3.11) equals 

the sign of the second derivative g" (z ) . 

sp 

b. Numerical Integration Using Newton-Cote's Formulas 
This standard technique [Ref. 11] is the method 
used by the program to perform a numerical integration over 
a subinterval of the total region of integration. Subinter- 
vals can be defined by dividing the total region of integra- 
tion into equal parts or by a scheme as described in the next 
section. Each subinterval is integrated by evaluating the 
Simpson's rule estimate: 



S = •=■ ( s + 2s. + 4 s,. + 2 s + 4s. + 

jo 1 2 j 4 



. + 2s . + s ) 
n-1 n 



(3.13) 
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where s 



denote values of the integrand evaluated at 
different points of the subinterval, and the distance between 
these points is given by the step size h. The original step 
size h is equal to half the size of the subinterval to be 
integrated. The step size is then halved to obtain a more 
accurate estimate which takes more points into account. This 
is repeated until a tolerance criterion is met. Using succes- 
sive Simpson's rule estimates, one can do a Richardson's 
extrapolation [Ref. 12]. This extrapolation uses the last 
two estimates by Simpson's rule, S£ and S^, where S 2 uses 

half the step size of S^. Both estimates have a global error 
4 

of order h . The extrapolated value becomes: 

Extrapolated value = S 2 + ( S 2 -S^)/15 (3.14) 

with an error of order h^. If two successive extrapolated 
values are within tolerance, the last extrapolated value is 
taken as the value of the integral for that subinterval. This 
scheme creates a semi-adaptive integration procedure because 
different subintervals can require different step sizes. It 
could be further enhanced but would then require more overhead. 
Another possibility would be the use of a Gaussian quadrature 
rule which needs less function evaluations for the same order 
of error but has the disadvantage that the points are not 
identically spaced and the scheme cannot be made adaptive in 
a computation efficient manner. 
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c. Numerical Integration Using the Euler Summation 

A problem in evaluating an integral like equation 
(3.4) is convergence. The real and imaginary parts of the 
integrand of equation (3.4) are oscillatory. Successive 
positive and negative contributions tend to cancel, but each 
contribution itself is not negligible and the integral behaves 
oscillatory as a function of its limits of integration. There- 
fore, one has to integrate over a relatively large region to 
obtain convergence. This is in conflict with the desired 
computational efficiency. 

Squire [Ref. 13: pp. 172-173] describes an inte- 
gration procedure for oscillatory integrands which splits the 
range of integration into parts using the zero crossings of 
the integrand as points of separation. The individual sub in- 
tervals can be evaluated with the. use of standard methods 
(as described in III.A.l.b) and summed to form the value of 
the integral. A technique like the Euler summation [Ref. 13: 
pp. 172-173] can be used to improve convergence of this sum- 
mation. The Euler summation obtains from an original sequence 




n 



a second sequence 



b 



1 



a 2 + a-^ , b 2 




b k * a k +a k +1 



and a third sequence 
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C 1 = b 2 + b l' c 2 = b 3 +b 2 c k = b k + l +b k 

and so on. It can be shown that the sum of 

2 a l + 4 b l + 8 C 1 + "• + ^in m l 

is the same as the sum of the original sequence. Convergence 
however is much improved which cuts down on the number of 
subintervals to be integrated and improves efficiency. 

This method works well on integrals where successive 
contributions from subintervals to the integral are of decreas- 
ing magnitude and of alternating sign as illustrated in Figure 
2. To use this method it must be possible to solve for the 
zero crossings of the integrand, which generally requires a 
closed form expression to solve for the roots of the integrand. 

2 . Application of Integration Methods 

The methods introduced in the previous section are 
used to define a procedure to calculate H(f,i,q). First, 

the integration w.r.t. direction cosine u is considered, 

J o 

which is given by 

/ S (f,u )exp{-jk [u AX . + (1-u^-v^ ) b//, ^AZ] }du (3.4) 

x o ^ J n o l o o o 

Both the method of stationary phase and the numerical 
method (Sections III.A.l.b and III.A.l.c) using the Euler 
summation are implemented. To apply the method of stationary 
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INTEGRAND 
G ( Z ) 




Fig. 2 Oscillatory Integrand Showing Decreasing 
Contributions to the Integral 
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phase we define in accordance with equation (3.8) 



g (u ) = -k [AX . u + (1 -u 2 - v 2 ) 1/2 AZ] (3.15) 

J o n i o o o 

The amplitude term S (f,u ) is the directivity function which 

x o 

is indeed slowly varying w.r.t. g(u Q ) • We find the sta- 
tionary point UOSP by setting g'(UOSP) equal to zero: 

g'(UOSP) = -k AX. +k AZ(1 -u 2 - v 2 )" 1/2 u = 0 (3.16) 

J n l n o o o 

leading to 

AX. 1 -v 2 1/2 

UOSP = -r^( ^ 2 ) (3.17) 

A 1 + (AX../AZ) z 

For the values of the second and third derivatives one finds 



g" (u ) = kAZ (1 -u 2 - v 2 ) 3/2 (l -v 2 ) 

o n o o o 

g"' (u ) = 3k AZu (1 - u 2 -v 2 ) -5//2 (l - v 2 ) 
o n o o o o 



(3.18) 



(3.19) 



We now form the ratio of the third and second derivatives 
at u q equal to the stationary point UOSP: 



g"' (UOSP) 
g" (UOSP) 



2 /2 (AX 2 + AZ 2 ) 1/2 

3AX. (l-vj 1/Z ^ 

1 ° AZ 2 



(3.20) 



The condition in equation (3.12) becomes 
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<< 1 



(3.21) 



AX. 

5 (_ AZ' ) 



>1 + ( AX . /A Z ) 2 



1 - v 



from which we conclude that the above condition is obeyed 

only when AZ >> AX^, which means that the method is only 

valid for relatively small cross ranges x^_ - x q since 

AX . = x -x +i d 1 . 
l r o x 

The program implements the following procedure to 
compute the magnitude and phase of the integral w.r.t. u q 
[ see equation (3.4)], denoted by IWRTUO. 

1) Calculate the stationary point UOSP 



UOSP = 



AX. 

l 

AZ 



1 - v' 



1 + (AX^/AZ) 



(3.17) 



2) Calculate the value of the second derivative DRV2ND 
at the stationary point: 



DRV2ND = 



k AZ (1 - UOSP 2 - v 2 ) 3/2 

o o 



d -vj) 



(3.22) 



3) Calculate magnitude and phase of the integral: 



MAG = S x (f,UOSP) /2 tt/DRV2ND 



(3.23) 



PHASE = -k [AX. UOSP + (1 - UOSP 2 - v 2 ) 1//2 AZ] +x 
n l o 4 



Note that the integral is a function of v q and the indices i 
and n through k n and AX^ . The dependence of the magnitude 
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on the index n is negligible, thus DRV2ND is calculated for 
all n using k Q . We will denote the integral by IWRTUO ( v q , k n , i) . 

The numerical method using the Euler summation can 
also be applied to the integration w.r.t. u . This numerical 
method is implemented to verify the result of the stationary 
phase method and extend application of the program to larger 
cross ranges. The applicability of the method can be shown 
by separating the integrand of equation (3.4) into a real 
and imaginary part 



/ (S ( f , u ) cos { -k [u iX. + (1 - u^ - vV^AZ]} 
x o n o i o o 

+ jS (f,u )sin{-k [u AX. + (l-u 2 -v 2 )''" //2 AZ]})du 



(3.25) 



n o l 



o o 



As shown above, the point of main contribution to the 
integral in equation (3.25) is the stationary point UOSP. 

Away from this point the integrand becomes more and more 
oscillatory and consequently more difficult to integrate 
numerically. One would therefore like to use a scheme in 
which the size of the subintervals is dependent on the behavior 
of the integrand. The method described in III.B.l.c uses 
the zero crossings of the integrand as points of separation, 
thus making sure that the integrand within a subinterval is 
"well-behaved" and can be efficiently integrated by standard 
numerical methods. This division into "well-behaved" subin- 
tervals could also be obtained by using the relative extremes 
as points of separation. The integrand of equation (3.25) 
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has closed form solutions to the zero crossings and relative 
extremes of both the real and imaginary parts. The zero 
crossings of the imaginary part of the integrand are found 
by solving: 



of the cosine term appearing in the real part of the integrand 
In order to avoid integrating the real and imaginary parts 
of the integral separately, the same subintervals are used 
for both the real and imaginary parts of the integral. Points 
of separation for the subintervals are the zero crossings of 
the imaginary part of the integrand, which coincide with the 
relative extremes for the real part of the integrand. As a 
starting point for the integration, define the first subinter- 
val using the point of main contribution, i.e., the stationary 
point UOSP . 



ary parts of the integrand around the stationary point. The 
figure also indicates the position of the subintervals. 
Integration of the successive subintervals (from UOSP outward) 
will form alternating terms of decreasing magnitude for both 



-k AX.u - k AZ (1 -u 2 -v 2 ) 1/2 
n i o n o o 



niT ; n integer (3.26) 



or 



u 



o 




(3.27) 



These also are the locations of the relative extremes 



Fig. 3 shows a typical behavior of the real and imagin 
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phase 



Fig. 3 Typical Behavior of Integrand in 

Eq. (3.25) Around Stationary Point 
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the real and imaginary parts of the integral, making it 
possible to effectively speed up the convergence by using a 
complex Euler summation. The decreasing magnitude is caused 
by the non-linear behavior of the phase term and the decreas- 
ing amplitude of the integrand which is given by the transmit 
far-field directivity function S x (f,u Q ). The program imple- 
ments the following procedure: 

1) Calculate the value of the stationary point UOSP 
according to equation (3.17). This value of u also 
denotes the maximum (negative) value of the phase 
term. 

2) Calculate NZ the maximum integer multiples of tt 
contained in the maximum (negative) value of the 
phase term (see Fig. 3) . 



NZ = INTEGER ( 



-k [UOSPAX. + (1-UOSP 2 -v 2 ) 1//2 AZJ 
n l o 



) (3.28) 



Note that NZ denotes the zero crossings of the imaginary 
part of the integrand nearest to the stationary point 
UOSP. 

3) Calculate the values of u Q for the zero crossings of 
the imaginary part of the integrand denoted by NZ 
using equation (3.27). 

AX . (NZtt) ±AZ [ (l-v 2 )k 2 (AX 2 +AZ 2 ) - (NZtt) 2 ] 1/72 
i o n i 

k (AX 2 +AZ 2 ) 
n i 

(3.29) 

Integrate the subinterval between these two values 
of UOZ. 

4) Decrement NZ and calculate the values of UOZy 2 > 
denoting the next zero crossings of the imaginary 
part of the integrand, .by using the formula given 
in 3) above. 



UOZ 



1,2 
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5 ) 



Define two subintervals, one on each side of the 
stationary point, from zero crossing to zero crossing 
for the imaginary part and from extreme to extreme 
for the real part of the integrand. 

6) Integrate both subintervals and add the contribution 
to the total value of the integral using a complex 
Euler summation. 

7) Check for convergence by comparing the difference 
between the last two values of the Euler summation 
with a tolerance. If not convergent, go back to 4) 
and repeat the procedure. 

One more note can be made concerning the numerically 
calculated value of the integral w.r.t. u for different k . 
Recalculation of the whole integral for different values of 
n would be very time consuming. We therefore introduce a 
correction term to obtain values of the integral for n ^ 0 , 
based on the value of integral for n = 0 . This correction 
term is only of importance to the phase of the integral, the 
effect of different values for k n on the magnitude is 
neglected. For n = 0 the stationary phase method results in 
the following expression for the phase: 



PHASE (n = 0) = -k [AX.UOSP+(1-UOSP 2 -Vo) 1//2 AZ]+ £ (3.30) 

O 1 u 4 



Brekhovskikh [Ref. 14] demonstrates that taking into 
account higher-order terms in the method of stationary phase 
does change the amplitude but not the phase of the resulting 
approximation. This shows that the point of main contribu- 
tion, i.e., the stationary point determines the phase of the 
integral. If the method of stationary phase is no longer 
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valid we can write a general form for the phase based on 
equation (3.30) : 

PHASE (n = 0 ) = -k Q [AX i UOSP+(l-UOSP 2 -v^) 1//2 AZJ + j + A<J> 

(3.31) 

where A<p represents some unknown correction term. This relies 
on the fact that the stationary point UOSP is always very 
close to the point of main contribution. Based on this, the 
phase term for n / 0 can be written as: 

PHASE (n) = - [k +(k -k ) ] [AX . UOSP+ (l-UOSP 2 -v 2 ) 1//2 AZ] +x + At}) (3.32) 
o n o l o 4 r 

where 



(k -k ) [AX. UOSP + (l-UOSP 2 -v 2 ) 1</2 AZ] 
n o i o 



(3.33) 



represents the correction term which is readily computed for 
different values of n. This phase correction term has to 
be added to the phase for n = 0 to obtain phase values for 
n / 0 . Note that this correction term is only important for 
larger ranges AZ, assuming that AZ >> AX. This can be seen 
by calculating a typical value for k n -k Q , where the element 
spacing is taken as A/2 and the sound speed profile has a 
constant gradient g: 

k -k = 2iTf( — — ) - - 12 n = 3.6 x lo ' 5 n (3.34) 

no c c c 
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where we used c q = 1500 m/sec and g = 0.017 1/sec. An upper 

limit for the correction term is (k -k )AZ if AZ >> AX. 

n o 

Having established two methods to implement the calcu- 
lation of the inner integral w.r.t. direction cosine u , the 
outer integral w.r.t. direction cosine v q is now considered. 

The integral can be written as [see equation (3.6)]: 

1 (N-l)/2 

I d n H M (k n' v o' n ' q)expl ‘ jk n v o aY qn 1IWRTU0(v o' k n' i)dv o 

a n=- (N-i ) / Z ^ 

q 

(3.35) 

The ocean medium transfer function values are calculated with 
the use of equations (2.11) , (2.16) and (2.18) . Note that the 

deterministic and random phase terms pertain to the case of 
a sound speed profile with a constant gradient. In applying 
these formulae, we have to take into account that the source 
depth for a transmit element is not y Q but y Q +ndy . This results 
in the following formulae which were actually implemented: 



A.. = (k v ) 

M o o 



- 1/2 



(3.36) 



n 



k c 

0 MD (n,q) ” 2 v [ g c (y +qd 

o ^ J r ^ y 



rr - 1) + iY qn I 



(3.37: 



k c 
n r n 



e MR (n ' < 3> - (f [ T on NR lnl c(y +qd'> n 

o ^ - 1 r ^ y 



(3.38! 



We note that both the amplitude and phase terms 



depend on the value of the wavenumber k n and direction cosine 
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v q . The expressions within the square brackets for 0 MD and 
0 ^ are only dependent on the geometry and can be precomputed. 
A random number generator was used to generate a series of 
N(0,1) distributed numbers to be used for n NR in equation 
(3.37) . 

Because of the complexity of the total phase term in 

the integrand of equation (3.35), a closed form approximation 

could not be found. The numerical method used is also based 

on the procedure given in III.A.l.c. The integral is divided 

into subintegrals in such a way that the subintegrals are 

all "well behaved." The problem is to find a general form 

which approaches the phase of the integrand. This form can 

then be used to calculate the points of separation in order 

to define the subintervals. We will proceed by separating 

the phase term of the integrand into a major part and a minor 

slowly varying part. Combining the phase term of the integral 

w.r.t. u q given by equation (3.32) with equation (3.35) 

results in the following expression for the integration w.r.t. 

v : 
o 



1 (N-l)/2 

a n=-(N-l)/2 nMno noy y4 



exp{-j (k -k ) [AYv +AX.UOSP+(1-UOSP 2 -v 2 ) 1//2 AZ] } * 
n o o i o 



IWRTUO(v ,k ,i) |exp{-jk [v AY+UOSPAX-+(l-UOSP 2 -v 2 ) AZ] }dv (3.39) 

on'^^oo i o o 
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where Ay = y - y . If we let 
1 r 1 o 



(N-D/2 

X(vJ = l d n H M (k n/ v o ,n,q)exp{-jk n v o (qd' -nd ) +j+ 



n=- (N-l)/2 



exp{-j (k n -k Q ) [ AYv_+AX^UOSP+ ( l-UOSP 2 -v 2 ) 1//2 AZ] } • | IWRTUO(v o ,k n ,i) | (3.40) 



and 



2 2 1/2 

a ( v ) = -k [AYv +UOSPAX. +(l-UOSP -v ) ' A Z ] 

o o o 1 o 

(3.41) 

a ( v ) = -k [AYv +(l-v 2 ) 1/2 Ax 2 + A Z 2 ] 

0 o o o 1 

Then the integral given by equation (3.39) becomes: 

1 

/ X ( v q ) exp [ j a ( v ) ] dv Q 

ci 

q 

where X(v q ) is slowly varying compared to exp{ja(v Q )}. The 
phase term a(v ) can now be used to define subintegrals which 
are "well behaved." This is accomplished by taking the roots 
of sin(a(v Q )) or cos (a(v Q ) ) as the points of separation between 
the different subintervals. Contributions from the subintegrals 
are in general of decreasing magnitude and alternating sign, 
so the Euler summation can be used to effectively speed up 
convergence . 

The implemented procedure uses the roots of sin(a(v Q ) ) 
as the points of separation and is similar to the one already 
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The point of main 



described for the integration w.r.t. u . 
contribution is approximately where the phase term a (v ) is 
stationary. This point is used to begin the numerical inte- 
gration. The stationary point VOSP is given by 

VOSP = AY/(R^ i + AY 2 ) 1//2 (3.43) 

which is a solution of a' (VOSP) = 0 where R . is the radial 

n 

distance from the center of the transmit array to a receiver 
element with index i and is given by 

2 2 1//2 

R . = (AX + AZ ) (3.44) 

ri l 

The points of separation for the integral can be found by 

1/2 

-AY(mr)±R . [k 2 (AY 2 +R 2 . ) - (nu) 2 ] 

VOZ = — 5 (3.45) 

k o (R ri + AY } 

Another method using a brute force technique was used 
to validate the above described procedure. This method 
integrates w.r.t. direction cosine v q between user specified 
limits of integration. It divides the range of integration 
into subintervals of equal size which are integrated by a 
standard routine as described above in III.A.l.b. Both 
methods agreed numerically although the method using the 
unequally spaced subintervals and the Euler summation was 
much more efficient in terms of computation time. 
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B. PROGRAM OUTLINE 



The flowchart in Fig. 4 gives an overview of the program 
logic. This logic is controlled by flags which are read 
from the input data file together with all program parameters. 
The program parameters determine: 

- Model geometry: x , y , z , x , y , z , M, N, M', N', 

d , d , d ' , d ' . ° ° ° r r r 

x' y x y 

Speed of sound profiles: c , g, a. 

Beam steering angles for the transmit array and the 
frequency to be used for the beam steering calculations. 

The fundamental period T 0 of the electrical input signal 
and the total number of samples L in the electrical 
output signal. The program determines the sampling 
period T s = T 0 /L and calculates the output electrical 
signal at time instants 



-(L-1)T q / 2, ..., 0, ..., (L-1)T q /2 



The values of the frequency components in the electrical 
input signal and their complex Fourier series coefficients 

c k* 

Tolerance values ERRORU and ERRORV for the numerical 
integrations w.r.t. u and v . The user should validate 
his choice for the tolerance values by comparing them 
against the absolute values calculated by the program 
for the integrations. 

Processing is done in accordance with flag settings and 

allows for a number of options : 

The frequency response H(f,i,q) can be calculated or 
read from a previously created file. 

The calculated frequency response H(f,i,q) can be 
filed for later use and/or printed. 

- The output electrical signal can be calculated and 
filed for further processing by other programs. 
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Fig. 4 Flowchart of Main Program Module 
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Most important numerical calculations required by the 
model are done within the module CALCHF. The module CALCHF 
calculates the frequency response matrix H(f,i,q), i.e., the 
frequency response for all receiver elements (i,q). This 
mainly involves evaluation of the double integral w.r.t. u 

o 

and v . It has two versions dependent on the integration 
technique used to evaluate the integral w.r.t. u q : 

CALCHF_1 Use stationary phase method to integrate 

w.r.t. u . 

o 

CALCHF_2 Use numerical technique to integrate w.r.t. u q . 

To avoid recalculation and improve program efficiency 

the outer integration w.r.t. v is done simultaneously for 

o 

all elements with the same index i. This is particularly 
important when using the more time consuming numerical method 
to integrate w.r.t. u q . This can be done by separating the 
summation inside the integrand into a part dependent on q 
and a part dependent on i. Recall that equation (3.35) 
gives the form of the integrand as follows: 

(N-l)/2 

d n H M (k n ,V o ,n,q)e ^ {_jk n V o AY qn' I ' ,JRrU0(v o ,k n ,i) (3.46) 

So given a value for IWRTUO for a given index i we can evalu- 
ate the integrand for all values of q, thus avoiding recalcu- 
lating the most computation intensive part. This scheme does 
result in more complex coding. The program has to do N' 



l 

n=- (N-l) /2 
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integrations and keep track of N' Euler summations simul- 
taneously. It does, however, also reduce some overhead. 
Convergence checking on the integration w.r.t. v q is done 
only on the two outermost elements with indices q = (N'-l)/2 
and q = - (N ' -1) /2 . 

The calculation of those terms that multiply IWRTUO(k ,v ,i) 
in the summation given by equation (3.46) is simplified by 
precomputing part of the phase components of the ocean medium 
transfer function given by equations (3.36) and (3.37). The 
program does this in the subroutine HMWKB and stores the 
resulting phase integral values in an array. To obtain a 
particular value for the phase, the program only needs to 
multiply with the factor k n /v Q . Especially in the case of a 
more complicated sound speed profile, the computation of 
the phase integrals can be involved and precomputation im- 
proves efficiency. The stated time-invariance of the model 
also forces precomputation since the random numbers used to 
model the phase term 0^ can be drawn only once. The random 
numbers were obtained from the IMSL routine GNNML which was 
used to generate a series of N xN 1 random numbers with distri- 
bution N(0,1). From these random numbers the different random 
phase terms 9 for all combinations (n,q) are readily com- 
puted using equation (3.37). Figure 5 shows a simple flow- 
chart of the CALCHF module. 

The innermost block which actually calculates H(f,i,q) 
is shown in more detail in Fig. 6. The diagram shows how the 
subintervals are defined by the program. The procedure followed 
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Fig. 5 Flowchart of Module CALCHF 
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V0SP: Approximate point of main contribution of the integral. 
NZ : Maximum negative integer denoting the multiples of 

contained in the maximum negative value of the phase. 
Note: The "Integrate Subintervals" block makes a Subroutine 
call to evaluate the integrand given by Eq. (3.46). 



Fig . 6 



Calculation of Subintervals by CALCHF 



is described in III. A. 2 for the integration w.r.t. v q and is 

similar to the procedure used for the numerical integration 

w.r.t. u in Section III. A. 2. 
o 

Implementation of the integration over a subinterval using 
the Newton Cote's formulas is straightforward and will not 
be examined in detail. The framework for this part of the 
program can be found in Squire [Ref. 13: p. 281]. 

The Euler summation is used to speed up convergence of 
the numerical integration. The objective is to form the sum 




from the original sequence ..., a n which are the 

contributions from the subintervals 1 through n, where 

b k = a k + a k+l ' c k = b k + b k+l ' 

The algorithm used is described below and uses an array 

COEF to store the coefficients a,b . , c _ , d etc. 

n n-1 n-2 n-3 

The size of this array determines how many contributions can 
be handled. 

DIMENSION COEF (2, SIZE) 

REAL EULSUM, FACTOR 
INTEGER NEW, OLD 
C Add contribution N 

COEF ( NEW , 1 ) = N th CONTRIBUTION 
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DO 100 J = 1, N 

100 COEF (NEW, J+l) = COEF (NEW, J) + COEF(OLD,J) 

FACTOR = FACTOR / 2 

EULSUM = EULSUM + COEF (NEW , N) /FACTOR 
C Prepare for next contribution 

SWAP ( NEW, OLD ) 

The program implements this algoritm in complex form to 
perform N' summations simultaneously. 

More details on the actual programming of the discussed 
numerical methods and flowcharts, data structures used, etc., 
is considered applied programming and will not be discussed 
in this thesis. 
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IV. COMPUTER SIMULATED DATA 



This section describes the validation of the model. It 
presents and interprets the data generated by the computer 
model for a number of test cases. 

A. PRESENTATION OF COMPUTER SIMULATED DATA 

Data generated by the computer model consists of the 
frequency response matrix H(f,i,q) and the output electrical 
signal at each element y(t,i,q) . Interpretation of the data 
is done directly by examining the frequency response H(f,i,q) 
and by the use of a three-dimensional DFT program to calculate 
the Fourier transforms of the computer simulated output elec- 
trical signals w.r.t. time and space. The program is based 
upon 3-dimensional DFT beamforming for planar arrays as des- 
cribed by Ziomek [Ref. 4]. First a brief overview of the 
output from the 3-D DFT program is given. 

The program calculates the 3-dimensional DFT of the 
electrical output signal from all receiver elements and is 
given by 



M' N' L' 



Yg (q,r,s) 



I l l 



l c m d n y (X.,m, n )exp[-j2T r qV L ] * 



m=-M' n=-N ' £=-L' 



exp[j2TTrm/(M-fMz) ]exp [j2irsn/(N+NZ) ] 



(4.1) 
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where : 



L: total number of time samples in one fundamental 

period T q of the signal. 

L' = (L-l)/2 

T g : sampling period 

M: total odd number of elements in the x-direction 

in the receive array 

M' = (M-l)/2 



N: total odd number of elements in the y-direction 

in the receive array 



N’ = (N-l)/2 

d , d : spacing in x- and y-directions , respectively, 

x y 

c , d : complex separable weights for the array 

given by 



c 

m 



a exp{j0 
m c J 



m 



} 



d 



n 



t> n e xp{ 



q,r,s: binnumbers of the DFT which are related to 

values of frequency f, direction cosine u 
and direction cosine v, respectively. 



y(£,m,n): output electrical signal at time instant 

£T s at element (m,n) of the receive array. 

MZ , NZ : Number of zeros used for padding to increase 

the resolution of the direction cosines u 
and v, respectively. 



Note that the symbols used to characterize this planar receive 
array are not the same symbols used to describe the receive 
array in the computer model (II. A. 1). Rectangular amplitude 
weighting without beam steering is used for all the results 
presented in this chapter. 



54 



In II. B. 3 we made the constraint to represent the transmit 
electrical signal by a finite Fourier series with maximum 
frequency KMAX • f . Thus, to avoid aliasing, we must sample the 
received signal with a sample rate equal to or greater than 
the Nyquist rate: 

f > 2 KMAX f T >2 KMAX T (4.2) 

s — o o — s 

The total number of samples L to be taken in order to avoid 
aliasing must obey the following inequality: 

L >_ 2 KMAX + 1 (4.3) 

The sampling period is determined by: 

T = T /L (4.4) 

s o 

The program evaluates the 3-dimensional summation in 
three steps. The first step is to perform a DFT w.r.t. time 
at each receiver element (m,n) . If the sound field incident 
upon the array is a general plane wave, then the received 
electrical signal can be expressed by 

y ( £ ,m,n) = g ( £T - t ) 

so 

where t = (u md + v nd ) /c denotes the time delay between 
o o x o y 

the signals at the different elements. The DFT w.r.t. time 
becomes 
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Y g (q,m,n) 



(4.5) 



c d 
m n 



Z=-L 



a {IT -t 
v s o 




where is the factor exp{- j 2Trq£/L} . Using the above stated 

J-i 

fact that the received signal is given by a finite Fourier 
series, we obtain 



Yg (q,m,n) = c m d n Lc q ex P {- j 2irqf o ( U Q md x + v Q nd y ) /c } (4.6) 

A normalized expression is defined by 

Yg N (q,m,n) = Y g ( q ,m , n , ) / ( C m d n L ) 

= c exp{-j2irqf (u md + v nd ) /c } (4.7) 

q^-^^ooxoy 

So Yg N (q,m,n) allows us to estimate the amplitude of the fre- 
quency components in the signal at each element (m,n) . 

The next DFT calculation is: 



m ' 

Yg ( q , r , n) = £ Y ( q , m , n) exp{ j 27rmr/ (m+mZ ) } (4.8) 

m=-M ' 



where we have padded with NZ zeros to increase resolution. 
Substitution equation (4.6) into equation (4.8) leads to 



Y~(q,r,n) = Led exp{-j2TT qf v n d /c } * 



q n 



o o y 



M' 



,rm 



J c exp{-i2Trqf u md /c}W., 

L ., , m ^ J 1 oo x / M+MZ 
m=-M 



(4.9) 
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and a normalized expression can be written as 



Y SN (q,r,n) = exp{-j 2tt q n d„/c} * 



o o y 



M' 



M' 



l c m exp{-j2Tr qf^ u^m d v /c}W M4 . M7 / \ 



m=-M ' 



'o o x' M+MZ / L . m 

m=-M 



(4.10) 



This expression is directly proportional to the Fourier 
series coefficients c and to the normalized far-field beam 

q 

pattern at a frequency q-f Q of a line array consisting of M 

elements which is lying along the X-axis. The planar array 

has N such line arrays as indicated by the index n. If no 

beam steering is done (c = a ) , then Y„..(q,r,n) has a maximum 

approaching c^ at the bin number related to the direction 

cosine value which is closest to u . The DFT bin number r 

o 

is related to the direction cosine value by 

u = rA/[(M+MZ)d ] ; A = c/qf (4.11) 

X o 



This value of u is an estimate of u . Therefore, for each 

frequency component of a plane wave propagating in the u q , v q 

direction, the expression given by equation (4.10) can be used 

to estimate the magnitude of the frequency component and the 

direction cosine u . An estimate of v.can be obtained in the 

o o 

same way by evaluating 
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Y SN (q,m,s) = c q ex P { ~ j 2^ q f Q u 0 m d x // ' c ^ ‘ 

N' N' 

I a expH 2,,f » nd /ctf ; / I b < 4 - 12 > 

n=-N' J n=-N ' 

The above mentioned estimates of c , u and v are provided 

q o o ^ 

by the program in the form of printouts of Yg^(q,m,n), 
Yg N (q,r,n) and Yg N (q,m,s). The printouts are done for user 
specified values of m and n. No padding with zeros is done 
for the transform w.r.t. time. Padding for the spatial trans- 
forms is done to obtain a user specified step size in direction 
cosine u or v. Graphical output for estimation of u q and v q 
was also obtained. The complete three-dimensional transform 
Yg(q,r,s) was used to produce a three-dimensional plot to aid 
in interpretation and for illustrative purposes. The esti- 
mates of u and v can also be used to calculate estimates 
o o 

of the spherical angles denoting the direction of the sound 
source as measured from the center of the receive array. 

These angles are also referred to as target bearing and eleva- 
tion angles. The estimates are given by 

\p = tan 1 (v /u ) (4.13) 

o o o 

e = sin -1 ( [u 2 + V 2 ] 1/2 ) (4.14) 

o o o 

/S A 

where i^ 0 and 0 q are estimates of the target bearing and ele- 
vation angles, respectively. Note that an error in the 
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estimation of either u Q or v D or both affects both the target 
bearing angle and target elevation angle estimates . 

B . TEST CASES 

Most test cases run by the model are kept simple to aid in 
interpretation. All cases are run first in a deterministic 
homogeneous medium so that results could be easily interpreted 
and compared to known methods and approximations. All test 
cases are checked to obey the criteria required by the trans- 
fer function. 

The values of the basic parameters such as source depth 
y Q , gradient g and reference speed of sound c q were chosen in 
order to model the wave propagation between a source located 
on the SOFAR channel axis and a deep receiver. At moderate 
latitudes from 60 °S to 60 °N the speed of sound on the SOFAR 
axis ranges from 1450-1485 m/sec in the Pacific and from 
1450-1500 m/sec in the Atlantic [Ref. 15]. The depth of the 
SOFAR axis is typically between 1000 and 1200 m at moderate 
latitudes [Ref. 15]. The gradient below the SOFAR axis is 
given by g = 0.017 1/sec [Ref. 19] . Values for the standard 
deviation a of the random component of the index of refraction 
range from 0.25 x 10 4 [Refs. 17,18] to 1.0 x 10 ^ [Ref. 16]. 
These values were used in choosing representative parameters 
for the following test cases: 

1 . Homogeneous Cases 

As a baseline test case we define the deterministic 
homogeneous case HMGl : 
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c = 1475 m/sec. 

g = 0: constant speed of sound, 

a = 0: deterministic medium. 

M = 11, N = 11: a 11 by 11 element transmit array with 

rectangular amplitude weighting in both the x 
and y directions. 



M' 


= 5, N 


1 = 5 


: a 5 by 5 element receive array. 


X 


= 0.0 


m 


x =0.0 


m 


o 






r 




y o 


= 1000 


. 0 m 


y = 2000.0 
1 r 


m 


z 


= 0.0 


m 


z = 1732.05 


m 


o 






r 




d 

X 


= d = 

y 


d' = 

X 


d' = 0.7375 m: 

y 


spacing equal to • 


e 

r 


= 30°, 


ifj 

v r 


90 ° : spherical 


angles denoting the line 



of sight between the center of the transmit array 
and the center of the receive array. 

u^ = 0, = 0.5: line of sight direction cosine values 

from the center of the transmit array to the center 
of the receive array. 

u b = 0, v b - 0.5: direction cosine values used to steer 

the main lobe of the transmit far-field beam pattern 
towards the center of the receive array. 

KMAX = 1: single frequency component. 

f = 1000 Hz, T = 0.001 sec. 
o o 

c^ = a-^ = 1.0: complex Fourier series coefficient equals 1. 

L = 3: 3 times samples per fundamental period T . 

Based on the definition of test case HMG1 we define 



the following cases: 



HMG2 : Same as HMG1 except 



x 



r 



u 



r 



0 



r 



200 m. 


z = 
r 


1720.464 m 


u b = 0.1, 


v = 
r 


ID 

O 

II 

rQ 

> 


30.66°, ip 


= 78. 


69° . 
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HMG3 : 



HMG4 : 



HMG5 : 



HMG 6 : 



HMG7 : 



1 IMG 8 : 



Same as HMGl except 

= 1500 m, z = 866.0246 m 

j r r 

Same as HMGl except 

y = 2500 m, z = 2598.074 m 

r r 

Same as HMGl except 

d - d = d' = d' = .36875 m 
x y x y 

f = 2000 Hz, T = 0.0005 sec. 

o o 

Same as HMGl except 

x = 1224.745 m z = 1224.745 m 

r r 

u r = Uj 3 = 0.6124, v = = 0 . 5 

0 = 52.24° , ip = 39.23° . 

r r r 

Same as HMGl except 

d = d = d' = d' = 0.184375 m 
x y x y 

f = 4000 Hz, T = 0.00025 sec. 

o o 

Same as HMG2 except 

KMAX = 5: five frequency components. 

f = 1000 Hz, T = 0.001 sec. 

o o 

f av = 5000 Hz. 
ma.x 

f^ = 3000 Hz: frequency used for steering the 

transmit beam pattern. 

d = d = d' = d' = 1475/2 f = 0.1475 m 
x y x y ' max 

c^ = a^ = 1.0; k = 1,2,..., 5: all complex Fourier 

series coefficients equal 1 . 

L = 11 : 11 time samples per fundamental period T . 



2 . Deterministic Inhomogeneous Cases 

The following deterministic inhomogeneous cases 
incorporate a sound speed profile with a constant gradient. 

The definitions are based on those used for the homogeneous cases . 
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INHMGl, INHMG2, and INHMG5 : Same as HMGl , HMG2 and HMG5 except 



g = 0.017 1/sec . 

INHMG8 : Same as HMG8 except 

g = 0.017 1/sec. 

3 . Random Inhomogeneous Cases 

The following cases with an index of refraction com- 
posed of a deterministic and a random component are defined: 

RINHMGl , RINHMG2: Same as INHMGl, INHMG2 except 

o = 1.0 x 10 -4 

RINHMG8 : Same as INHMG8 except 

a = 0.25 x 10~ 4 



C. RESULTS 

This section will present and discuss the results from 
the computer model runs on the test cases as defined in Section 
IV. B. Some cases were run using both methods of integration 
w.r.t. u . The two methods showed no significant differences. 

First, some trivial checks were made using the various 
homogeneous cases. Note, that all test cases have the same 
line of sight direction cosine v . For the homogeneous cases 
with a single frequency component of 1000 Hz the amplitude 
and phase terms of the ocean transfer function are equal and 
given by: 



*M 



(k v ) 
o o 



- 1/2 



MD 



= 0 



(4.15) 
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Since those cases have the same ocean transfer function 



values, the difference in the magnitude of the frequency 
response H(f,i,q) for the various cases should only depend 
on the total range given by 



and the transmit far-field beam pattern. Cases HMG1 , HMG3 

and HMG4 also have the same line of sight direction cosine 

value u , so the influence of the transmit far-field beam 
r 

pattern is the same for these cases and magnitude differences 
between the frequency responses can only depend on the total 
range and should behave in accordance with the expected 
spherical spreading for a homogeneous medium. Table IV. 1 
shows cases HMG1 , HMG3 and HMG4 together with their total 
range and the magnitude of the frequency response at the center 
element of the receive array, i.e., |H(f,0,0)|. We conclude 

that the magnitude indeed falls off as 1/r, only the total 
range is of importance. 

Table IV. 2 shows cases HMG1, HMG2 and HMG6 which differ 
only in cross range - x q and line of sight direction cosine 
value u^ but have the same total range. In these cases only 
the far-field beam pattern can cause differences in the magni- 
tude of the frequency response. The table shows that as the 
cross range increases the magnitude falls off, which is con- 
sistent with the fact that when a beam pattern is steered 
towards end-fire (to larger values of u^ and v^) the 



R 



'total 




(y r - y Q ) 2 + <z r - z o ) 2 (4.16) 
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directivity index decreases due to the increasing beamwidth 
of the main lobe . 

Table IV. 3 shows the difference in magnitude of the fre- 
quency response at the center of the receive array between 
HMGl , HMG5 , and HMG7 . Cases with the same geometry but a 
different value for the single frequency component. This 
difference in magnitude is due to the frequency dependent 
amplitude factor inherent in the WKB approximation, which is 
not exact for a homogeneous medium. 

Table IV. 4 compares the magnitudes of the frequency 
response between the homogeneous cases HMGl, HMG2 , and HMG5 
and their inhomogeneous counterparts INHMGl , INHMG2 and INHMG5. 
We see that the magnitudes for the inhomogeneous cases are 
slightly lower, an expected phenomena because of the channel- 
ing effect in the inhomogeneous medium. In a medium with a 
constant positive gradient, the acoustic rays are bent upwards 
and the sound intensity decreases with increasing depth. 

This causes a lower magnitude of the acoustic signal then 
expected by spherical spreading only. 

We will now present numerical and graphical output from 
the three dimensional DFT beamforming program. The frequency 
transform for the single frequency component cases simply 
reveals the magnitude of the transfer function multiplied by 
the amplitude of the input electrical signal and is not shown 
here. Figures 7 through 10 show the magnitude of the spatial 
DFT w.r.t. x, i.e., Yg N (q,r,n) versus direction cosine u, and 
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the magnitude of the spatial DFT w.r.t. y, i.e., Y (q,m,s) 
versus direction cosine u, for the cases HMG1, HMG2 , INHMGl 
and INHMG2 . These graphs allow us to estimate the direction 
of propagation of the incoming plane wave by finding the 
coordinates of the maxima on the graphs. Tables IV. 5 to 
IV. 8 provide numerical values for the above mentioned DFT 1 s 
w.r.t. x and y. They show values of Yg N (q,r,n) and Yg N (q,m,s) 
around the maxima and make it possible to form more accurate 
estimates . 

/N /\ 

From the estimated values u and v , we can calculate 

o o 

A 

estimates of target elevation angle 0 q and target bearing 

/N 

angle Table IV. 9 gives the results for the different 

A /S 

cases. The values of the estimators u and v were determined 

o o 

by performing a second order interpolation on the data given 
in tables IV. 5 through IV. 8. It is seen that the estimates 
of target angles for the homogeneous cases are correct and 
equal the line of sight angles from the receiver to the trans- 
mitter. The estimated target angles for the deterministic 
inhomogeneous cases deviate slightly from the line of sight 
angles. This is due to the positive sound speed gradient 
which causes upward bending of the acoustic rays so the direc- 
tion cosine v q , directly related to the vertical angle of 
propagation of the plane wave, decreases as the wave travels 
along. The net effect for these cases is relatively small 
but noticable. 

The effect of randomness, characterized by a random component 
in the index of refraction with a standard deviation of 
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-4 

a = 1.0 x 10 is examined next. Figure 11 shows estimates of 

u and v for the case RINHMG1 and Fig. 12 shows estimates of 
o o 3 

v for cases RINHMGl and RINHMG2 . As can be seen, the esti- 
o 

mate of u q is not influenced, which follows readily from the 

fact that the randomness is only present in the y direction. 

The estimate of v as shown in Fig. 12 for the case RINHMG2 

o 

shows the same random fluctuations as the one shown for case 
RINHMGl in Fig. 11, since both cases used the same seed for 
the random number generator, resulting in the same random 
number sequence. Figure 12 also shows RINHMGl with a different 
seed for the random number sequence. Although both are only 
two possible samples they clearly show the strong effect of 
the randomness on the beamformer output and the deterioration 
of the quality of the estimates (see Table IV. 9) . 

Two three-dimensional plots are shown in Fig. 13, one for 
the deterministic inhomogeneous case INHMG2 and one for the 
random inhomogeneous case RINHMG2 . These graphs illustrate 
the kind of output that can be obtained from the DFT beam- 
former, although it is difficult to accurately read off 
numerical values. 

Finally, the cases HMG8 , INHMG8 and RINGHM8 , which involve 
several frequency components, are considered. The output 
for these cases is given in numerical and graphical form. 

, Table IV. 10 compares the Fourier series coefficients of the 
output electrical signal with the Fourier series coefficients 
of the transmitted electrical signal. The widely varying values 
are caused by: 
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1) The beam pattern of the transmit array is frequency 
dependent. The differences between the transmitted frequency 
components are relatively large and cause large differences in 
amplitude of the Fourier series coefficients. The maximum 
values are observed at f = 3000 Hz, i.e., the frequency for 
which beam steering was done. 

2) The deterministic inhomogeneous medium: As already 

shown in other test cases the amplitude of the received signal 
is frequency dependent due to the WKB approximation to the 
wave equation in an inhomogeneous medium. Equation (2.10) 
clearly shows the frequency dependence of the medium transfer 
function . 

3) The random inhomogeneous behavior of the medium: The 

randomness introduces additional deviations. It can be seen 
from table IV. 10 that the effect of the randomness increases 
with increasing frequency. Also, in the case of a random 
medium, the estimation of the Fourier series coefficients 
becomes strongly dependent on the location of the element in 
the receive array. Note that table IV. 10 only shows values 
for the receive element (m = 0 , n = 0) . 

The received angular spectrum for the different frequency 
components of the cases INHMG8 and RINHMG8 is shown in Figs. 

14 through 19. For the deterministic inhomogeneous case, 

Figs. 14 through 16 show a regular behavior of the angular 
spectrum according to the rectangular amplitude weighting 
used for the receive array. Estimation of u q is exact. 
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Estimation of v q returns a slightly different value from the 

line of sight value for direction cosine v due to the 

inhomogeneous behavior of the medium. Figures 14, 15 and 

16 demonstrate the increasing directivity of the receive 

array with increasing frequency, thus leading to more accurate 

estimates for higher frequencies. 

Figures 17, 18 and 19 clearly show the strong influence 

of the randomness on the resulting angular spectrum. The 

randomness creates strong "false" sidelobes and causes an 

offset of the estimated value of direction cosine v . The 

o 

effects of the randomness increase with increasing frequency. 

At f = 5000 Hz (Fig. 19) the estimate of v is determined by 

o 

a spurious side lobe and returns a totally wrong value. 

Table IV. 11 shows estimated values of u and v at the differ- 

o o 

ent frequencies and their associated estimates of target 
elevation and bearing angles, 0 Q and ip Q , respectively. 
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TABLE IV. 1 



Magnitude of Frequency Response Comparing Cases 
Which Are Different Only in Total Range 



Case 


Total Range 
(in m) 


|H(f ,0,0) j 


HMG1 


2000 


0.02434 


HMG3 


1000 


0.04868 


HMG4 


3000 


0.01622 



TABLE IV. 2 

Magnitude of Frequency Response Comparing Cases 
Which Are Different Only in Cross Range 



Case Total Range Cross Range |H(f,0,0) 

(in m) (in m) 



HMG1 


2000 


o 

o 


0 .02434 


HMG2 


2000 


200.0 


0.02417 


HMG6 


2000 


1224 . 75 


0.01654 



TABLE IV. 3 



Magnitude 

Which 


of Frequency Response 
Are Different Only in 


Comparing Cases 
Frequency 


Case 


Frequency 
(in Hz) 


| H ( f ,0,0) | 


HMG1 


1000 


0.02434 


HMG5 


2000 


0.03441 


HMG7 


4000 


0.04866 



TABLE IV. 4 

Magnitude of Frequency Response Comparing 
Homogeneous and Inhomogeneous Cases 



Case 




I H ( f , 0 


,0) | 


HMGl - 


INHMGl 


0.02434 


- 0.02341 


HMG2 - 


INHMG2 


0.02417 


- 0.02326 


HMG5 - 


INHMG5 


0.03441 


- 0.03310 
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TABLE IV. 5 



Magnitude of 
Cosine 


Y SN (q,r,n) 
u for HMGl , 


Vs. Direction 
INHMGl 


Direction 
Cosine u 


HMGl 


INHMGl 


-0.015 


.0121416 


. 0116804 


-0 . 010 


.0121566 


.0116945 


-0.005 


.0121655 


.0117030 


0.0 


. 0121685 


. 0117058 


0.005 


.0121655 


. 0117030 


0.010 


. 0121566 


. 0116945 


0 . 015 


.0121416 


.0116804 





TABLE IV. 6 




Magnitude of 
Cosine 


Ysn ( q? r ,n) 
u for HKG2 , 


Vs. Direction 
INHMG2 


Direction 
Cosine u 


HMG2 


INHMG2 


0.085 


.0120593 


.0116000 


0.090 


.0120744 


.0116147 


0 . 095 


.0120836 


. 0116239 


0.100 


.0120869 


.0116274 


0.105 


.0120842 


.0116254 


0.110 


.0120756 


.0116178 


0.115 


.0120611 


.0116046 
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TABLE IV. 7 



Magnitude of 
Cosine 


ysn ( q,m, s) 

v for HMGl 


Vs. Direction 
, HMG2 


Direction 
Cosine v 


HMGl 


HMC-2 


0 .484 


.0121367 


.0120534 


0.489 


.0121531 


.0120696 


0.494 


.0121636 


.0120800 


0.499 


.0121681 


.0120844 


0.504 


.0121666 


.0120828 


0.509 


.0121591 


.0120754 


0.514 


.0121457 


. 0120620 





TABLE IV. 8 




Magnitude of 
Cosine 


Y SN (q,m,s) Vs. 
v for INHMGl , 


Direction 

INHMG2 


Direction 
Cosine v 


INHMG1 


INHMG2 


0.475 


.0116731 


.0115947 


0.480 


.0116893 


.0116111 


0.485 


.0117001 


.0116215 


0.490 


.0117051 


.0116264 


0 . 495 


.0117045 


.0116258 


0.500 


.0116983 


.0116196 


0 . 505 


.0116864 


.0126078 
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TABLE IV. 9 



Estimated Values u , v . 9 and ip 

o o o o 



Case 


u 

o 


/\ /\ 

v o 9 o 

(in degrees) 


A 

^ o 

(in degrees) 




HMG1 


0.000 


0.500 


+30.00 


+90.00 






HMG2 


0.100 


0 . 500 


+30.65 


+78. 69 






INHMG1 


0.000 


0.492 


+29 . 47 


+90 .00 






INHMG2 


0.101 


0.492 


+30.15 


+78.40 






RINHMG2 


0 . 101 


0.523 


+32 . 19 


+79 . 07 


( seed 


#1) 


RINHMG1 


0.000 


0.523 


+31.53 


+90.00 


(seed 


#1) 


RINHMGl 


0.000 


0.508 


+30 . 53 


+90.00 


( seed 


#2) 


1) The line of 


sight angles 


from receive 


to transmit 




array 


(or target angles) 


are : 








HMG1 , 


INHMG1 


, RINHMGl: 


0 r = +30.00°, 


\p = +90. 
r 


O 

O 

o 




HMG2 , 


INHMG2 


, RINHMG2 : 


0 = +30 . 65° , 

r 


\p = + 78. 
r r 


69° 





2) Note that equation (4.14) which is used to calculate 
estimates of target bearing angle from estimates of 
u Q and v Q returns two results. This is due to the 
ambiguity of a planar array. For the above table the 
appropriate values were selected. 
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TABLE IV. 10 



Estimation of Fourier Series Coefficients 
Y SN (q,0,0) for HMG8 , INHMG8 and RINHMG8 



Frequency 
(in Hz) 


HMG8 


INHMG8 


RINHMG8 


1000 


0.0010987 


0.0009705 


0.0010341 


2000 


0.0098126 


0 .0098901 


0.0102369 


3000 


0 .0209366 


0.0201906 


0.0177316 


4000 


0.0138747 


0.0122052 


0.0087726 


5000 


0 . 0024539 


0.0033176 


0 . 0051316 



Note: The Fourier series coefficients of the trans- 

mitted electrical signal are c. = a, = 1 ; 
k = 1,2, ... ,5 



TABLE IV. 11 



Estimated Values of u , v , 0 and ib 

o o o o 



Case 


Frequency 
(in Hz) 


A 

u 

o 


0 

< > 


A 

e o 

(in degrees) 


/\ 

(in degrees 


HMG8 


1000 


0 . 10 


0 . 50 


+ 31 


+ 79 


it 


2000 


0.10 


0 . 50 


+31 


+ 79 


it 


3000 


0.10 


0 . 50 


+ 31 


+ 79 


tl 


4000 


0.10 


0.50 


+31 


+ 79 


It 


5000 


0 . 10 


0 . 50 


+ 31 


+ 79 


INHMG8 


1000 


0.10 


0 . 49 


+30 


+ 78 


II 


2000 


0 . 10 


0 .49 


+30 


+ 78 


It 


3000 


0.10 


0 . 49 


+30 


+ 78 


It 


4000 


0.10 


0.49 


+ 30 


+ 78 


It 


5000 


0.10 


0 . 49 


+ 30 


+ 78 


RINHMG8 


1000 


0.10 


0.63 


+ 40 


+ 81 


It 


2000 


0.10 


0 . 49 


+ 30 


+ 78 


If 


3000 


0 . 10 


0 . 53 


+ 33 


+ 79 


It 


4000 


0.10 


0 . 52 


+ 32 


+ 79 


It 


5000 


0.10 


-0.43 


-26 


-77 



ESTIMATION OF UO 




ESTIMATION OF VO 




Fig. 7 Estimation of Direction of Propagation 
for Case HMGl 
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0.012 



ESTIMATION OF UO 




ESTIMATION OF VO 




Fig. 8 Estimation of Direction of 
Propagation for Case HMG2 
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MAGNITUDE MAGNITUDE 

O . Q 02 0.001 0.006 0.000 0.010 0.000 0.002 0.004 0.006 0.008 o.oto 



ESTIMATION OF UO 




ESTIMATION OF VO 




Fig. 9 Estimation of Direction of 
Propagation for Case INHMGl 
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ESTIMATION OF UO 




ESTIMATION OF VO 




Fig. 10 Estimation of Direction of 
Propagation for Case INHMG2 
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ESTIMATION OF UO 




ESTIMATION OF VO 




Fig. 11 Estimation of Direction of Propagation for 
Case RINHMG1 
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MAGNITUDE -10' MAGNITUDE *10 

0.0 l.o 2.0 3.0 4.0 5.0 6.0 7.0 0.0 0.0 t.O 3.0 3.0 4.0 S.O 8.0 7.0 8.0 



ESTIMATION OF VO 




CASE RINHMGl FREOUENCY : 1000 HZ 



ESTIMATION OF VO 




Fig. 12 Estimation of v for Cases RINHMG2 and 
RINHMG1 (Different Seeds) 
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Fig. 13 Angular Spectrum for Cases INHMG2 and RINHMG2 
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Fig. 14 Angular Spectrum for Case INHMG8 
at f = 1000 and f = 2000 Hz 




Fig. 15 Angular Spectrum for Case 
INHMG8 at f = 3000 Hz 
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Case RINHMG8 



W&ttYWOE. v .o 

U K 0 -& k 1 — 




Fig. 19 Angular Spectrum for Case 
RINHMGS at f = 5000 Hz 
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V. CONCLUSIONS AND RECOMMENDATIONS 



Computer simulation of a mathematical model of wave propa- 
gation through a random, space-variant, time-invariant, ocean 
medium between two planar arrays has been accomplished. The 
computer simulation incorporates an index of refraction which 
is only a function of depth. The deterministic component of 
the index of refraction represents a sound speed profile with 
constant gradient. The random component is assumed to have 
statistical properties independent of depth. The effects of 
the randomness on the electrical output signals at the differ- 
ent receiver elements are assumed to be uncorrelated. 

The computer program was written in FORTRAN. It implements 
both analytical and numerical integration techniques. The 
results from both techniques are in close agreement. 

A number of test cases have been run. Interpretation of 
the data for simple cases show results which are consistent 
with expectations from theory. The model correctly demonstrated 
the influences of the far field beam pattern of the transmit 
array, the source-receiver geometry, and the deterministic 
inhomogeneous behavior of the medium. Also, a random inhomo- 
geneous medium was shown to affect the wave propagation and 
distort the frequency and angular spectra, leading to errors 
in the estimation of the Fourier series coefficients, target 
bearing angle, and target elevation angle. A final test case 
RINHMG8, with a transmit signal containing several frequency 
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components demonstrated a complex dependence on the above 
mentioned factors, making interpretation of the results diffi- 
cult. The effect of the randomness was very pronounced in 
this case. 

The computer program could be further improved by imple- 
menting a more sophisticated method to simulate the randomness 
of the medium. The assumption that the random phase terms at 
the different receiver elements are uncorrelated breaks down 
for small interelement spacings. The computer program can 
also be easily modified to accommodate a more complex sound 
speed profile in which the gradient is a function of depth, 
for example. 

Future application of the model could be devoted to 
investigating the propagation of more complex transmit signals. 
One could also investigate the effects of the randomness of 
the ocean medium on signal processing using different sizes 
of receive arrays, i.e., increasing the number of elements. 



89 



LIST OF REFERENCES 



1. Ziomek, L.J., Underwater Acoustics — A Linear System Theory 
Approach , Academic Press, Orlando, Florida (unpublished). 

2. Ibid., Chapter 7, Section 7.2.1. 

3. Ibid., Chapter 7, Section 7.2.2. 

4. Ibid., Chapter 5, Section 5.1. 

5. Laval, R. , "Sound Propagation Effects on Signal Processing, 
Signal Processing , edited by Griffiths, P.L. and others, 
pp. 223-241, Academic Press, New York, 1973. 

6. Laval, R. , "Time-Frequency-Space Generalized Coherence 
and Scattering Functions," Aspects of Signal Processing , 
edited by Tacconi , G., Vol . I, pp. 69-87, D. Reidel 
Publishing Company, Dordrecht, The Netherlands, 1977. 

7. Laval, R. and Labasque , Y., "Medium Inhomogeneities and 

Instabilities: Effects on Spatial and Temporal Processing, 

Underwater Acoustics and Signal Processing , edited by 
Bjorno, L. , pp. 41-70, D. Reidel Publishing Company, 
Dordrecht, The Netherlands, 1981. 

8. Middleton, D. , "A Statistical Theory of Reverberation and 

Similar First-Order Scattered Fields. Part III: Waveforms 

and Fields," IEEE Trans. Information Theory, Vol. 18, 

pp. 35-67, 1972. 

9. Middleton, D., "The Underwater Medium as a Generalized 
Communication Channel , " Underwater Acoustics and Signal 
Process ing , edited by Bjorno, L. , pp. 589-612, D. Reidel 
Publishing Company, Dordrecht, The Netherlands, 1981. 

10. Officer, C.B., Introduction to the Theory of Sound Trans - 
mission , pp. 67-68, McGraw-Hill, New York, 1958. 

11. Gerald, C.T., Applied Numerical Analysis , pp. 60-72, 
Addison-Wesley , 1970. 

12. Moursund, D.G. and Duris, C.S., Elementary Theory and 
Application of Numerical Analysis, pp. 194-199, McGraw- 
Hill, 1967 . 

13. Squire, W., Integration for Engineers and Scientists , 
American Elsevier Publishing Company, New York, 1970. 



90 



14. Brekhovskikh, L.M., Waves in Layered Media , pp. 234-239, 
Academic Press, New York, 1980. 

15. Brekhovskikh, L. and Lysanov, Yu., Fundamentals of Ocean 
Acoustics , p. 3, Springer-Verlag , 1982. 

16. Clarke, R.H., "Sound Propagation in a Variable Ocean," 
Journal of Sound and Vibration , Vol. 34, p. 472, 1974. 

17. Mintzer, D., "Wave Propagation in a Randomly Inhomogeneous 
Medium I," Journal of the Acoustic Society of America , 

Vol. 25, p. 925, 1953. 

18. Mintzer, D., "Wave Propagation in a Randomly Inhomogeneous 
Medium I and II," Journal of the Acoustic Society of 
America, Vol. 26, p. 1110, 1953. 

19. Kinsler, L.E., and others. Fundamentals of Acoustics , 
p. 401, 3rd ed . , Wiley & Sons, 1981. 



91 



INITIAL DISTRIBUTION LIST 



NO 



1. Defense Technical Information Center 
Cameron Station 

Alexandria, Virginia 22314 

2. Library, Code 0142 
Naval Postgraduate School 
Monterey, California 93943 

3. Department Library, Code 61 
Department of Physics 
Naval Postgraduate School 
Monterey, California 93943 

4. Assistant Professor L.J. Ziomek, Code 62Zm 
Department of Electrical and 

Computer Engineering 
Naval Postgraduate School 
Monterey, California 93943 

5. Inspecteur Onderwijs Zeemacht 
Ministry of Defense (Navy) 

Koningin Marialaan 17 

Den Haag 

The Netherlands 

6. Vlagofficier belast met de of ficiersvorming 
Koninklijk instituut voor de Marine 

Het Nieuwe Diep 8 

Den Helder 

The Netherlands 

7. LT J. Vos 

p/a Dorpsstraat 12 
Schipluiden, 2636 CC 
The Netherlands 

8. Mr. Charles Stewart 
DARPA 

1400 Wilson Blvd. 

Arlington, VA 22200 



Copies 

2 

2 

1 

15 

2 

2 

3 

1 



9 2 



213187 



Thesis 

V9834 Vos 

c.l Linear time- invar- 

iant space-variant 
filters and the WKB 
approximation with 
applications to under- 
water acoustic signal 
processing. 







213187 


Thesis 

V9834 


Vos 


c.l 


Linear time-invar- 




iant space- variant 
filters and the WKB 
approximation with 
applications to under- 
water acoustic signal 
processing. 




